The runoutSIM package is developed for regional-scale
runout simulations (spatial distribution and velocity) using random
walks. This vignette provides an introduction to applying
runoutSIM for regional debris flow simulation using a
subcatchment of the upper Maipo river basin, Chile as an example. In
this vignette, we will:
Spatial handling in the runoutSIM package is supported
by the sf (vector data) and terra
(gridded/raster data) packages. In this example, we load/read a 12.5 m
spatial resolution DEM, debris flow runout polygons, and corresponding
source points mapped for each polygon. We also create a hillshade model
to help with visualization of the results.
Tip: Currently, this package is designed to work with data in a UTM coordinate reference system (CRS). You will need to project your own data to a local UTM CRS for it to work.
# load packages
library(runoutSim)
library(terra)
#> Warning: package 'terra' was built under R version 4.4.3
#> terra 1.8.54
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
# Load digital elevation model (DEM)
dem <- rast("C:/GitProjects/runoutSIM/Dev/Data/elev_fillsinks_WangLiu.tif")
# Compute hillshade for visualization
slope <- terrain(dem, "slope", unit="radians")
aspect <- terrain(dem, "aspect", unit="radians")
hill <- shade(slope, aspect, 40, 270)
# Load debris flow runout source points and polygons
source_points <- st_read("C:/GitProjects/runoutSIM/Dev/Data/debris_flow_source_points.shp")
#> Reading layer `debris_flow_source_points' from data source
#> `C:\GitProjects\runoutSim\Dev\Data\debris_flow_source_points.shp'
#> using driver `ESRI Shapefile'
#> Simple feature collection with 73 features and 1 field
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 389175.6 ymin: 6293926 xmax: 398788.6 ymax: 6327439
#> Projected CRS: WGS 84 / UTM zone 19S
runout_polygons <- st_read("C:/GitProjects/runoutSIM/Dev/Data/debris_flow_runout_polygons.shp")
#> Reading layer `debris_flow_runout_polygons' from data source
#> `C:\GitProjects\runoutSim\Dev\Data\debris_flow_runout_polygons.shp'
#> using driver `ESRI Shapefile'
#> Simple feature collection with 73 features and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 389139.1 ymin: 6293864 xmax: 398852.9 ymax: 6327455
#> Projected CRS: WGS 84 / UTM zone 19S
# Load basin boundary
bnd_catchment <- st_read("C:/GitProjects/runoutSIM/Dev/Data/basin_rio_olivares.shp")
#> Reading layer `basin_rio_olivares' from data source
#> `C:\GitProjects\runoutSim\Dev\Data\basin_rio_olivares.shp'
#> using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 15 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 386042.8 ymin: 6293348 xmax: 404405.3 ymax: 6331286
#> Projected CRS: WGS 84 / UTM zone 19SUsing the leafmap() function, we can quickly create an
interactive leaflet map for viewing the data. Similar to
dplyr, we use the %>% operator to pipe
together multiple leafmap() calls to create a combined
leaflet map with multiple layers.
Tip: leafmap() has the option to
interactively query raster values with a mouse hover. This option is
turned on by default, but can result in a considerable increase in the
file size of the leaflet widget (e.g. when exported for sharing). It is
recommended to have add_image_query = FALSE when saving /
exporting the leaflet map.
map <- leafmap(bnd_catchment, color = '#f7f9f9', fill_color = '#FF000000', weight = 4) %>%
leafmap(runout_polygons) %>%
leafmap(source_points, color = '#e74c3c') %>%
leafmap(hill, palette = grey(0:100/100), opacity = 1, add_legend = FALSE,
add_image_query = FALSE) %>%
leafmap(dem, palette = viridis::mako(10), opacity = 0.6, add_image_query = FALSE)
mapRunout paths are simulated from individual source cells, defined as
x, y coordinates. For example, if we are simulating runout from a single
source cell defined by a source point (vector data), we use
sf::st_coordinates() to extract the x, y values.
# Select a single debris flow and source point for the example
runout_polygon <- runout_polygons[31,]
# Get corresponding source point
source_point <- st_filter(st_as_sf(source_points), st_as_sf(runout_polygon))
# Get coordinates of source point
source_crds <- st_coordinates(source_point)
print(source_crds)
#> X Y
#> [1,] 395698.1 6307421Now that we have the coordinates, we can run the random walk runout
simulations using runoutSim(). Where:
mu and md are the sliding friction
coefficient and mass-to-drag parameters of the two-parameter friction
model controlling runout distance
slp_thresh, exp_div, and
per_fct are parameters controlling the dispersion of the
random walks
walks is the number of random walk
iterations
sim_paths = runoutSim(dem = dem,
xy = source_crds,
mu = 0.08,
md = 140,
slp_thresh = 35,
exp_div = 2.5,
per_fct = 1.95,
walks = 1000)
# plot structure of sim_paths
str(sim_paths)
#> List of 4
#> $ start_cell : num 2803625
#> $ cell_trav_freq: 'table' int [1:2012(1d)] 5 2 7 17 12 7 4 1 1 1 ...
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ all_indices: chr [1:2012] "2774166" "2774167" "2775633" "2775634" ...
#> $ cell_max_vel : num [1:2012] 15.8 16.2 12 14.3 16.2 ...
#> $ prob_connect : NULLrunoutSim() returns a list containing:
The cell number of the DEM indicating the location of the source cell
Traverse frequencies for each grid cell
Maximum velocities (m/s) for each grid cell
These values are stored as cell numbers (e.g.,
$all_indices). Although not calculated here,
prob_connect provides the probability the source cell
connected to a feature of interest.
To visualize the results, we convert the output from
runoutSim() to a raster.
## Convert paths to raster with cell transition frequencies
paths_raster <- walksToRaster(sim_paths, dem, method = "freq")
# Plot results
paths_raster <- crop(paths_raster, ext(runout_polygon)+500)
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
mar=c(2,2,1,4), main = "Runout traverse frequency")
plot(paths_raster, legend = T, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)In walksToRaster(), we use method = "freq"
to return traverse frequencies. We can also return traverse (ECDF)
probabilities by using method = "cdf_prob".
We use velocityToRaster() to get a raster of velocity
values (m/s).
## Convert paths to raster with max. velocities
vel_raster <- velocityToRaster(sim_paths, dem)
# Plot results
vel_raster <- crop(vel_raster, ext(runout_polygon)+500)
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
mar=c(2,2,1,4), main = "Runout velocity")
plot(vel_raster, col = map.pal('plasma'), legend = TRUE, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)The runoutSim() function also allows us to calculate the
connectivity probability of a source cell to feature of interest -
i.e. the proportion of walks that intersect a feature from a given
source cell.
In this example, we will determine the connectivity probability of a
simulated debris flow to a river channel. Using
makeConnFeature() function, we will create a matrix
indicating with a value of 1 which cells contain the river channel.
# Load river channel polygon (vector)
river_channel <- st_read("C:/GitProjects/runoutSIM/Dev/Data/river_channel.shp")
#> Reading layer `river_channel' from data source
#> `C:\GitProjects\runoutSim\Dev\Data\river_channel.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 393504.8 ymin: 6293330 xmax: 396806.5 ymax: 6328149
#> Projected CRS: WGS 84 / UTM zone 19S
# Create a connectivity feature for runoutSim that matches the input DEM
feature_mask <- makeConnFeature(river_channel, dem)Now that we have a connectivity feature, we can re-run
runoutSIM() with some additional parameter options:
source_connect = TRUE allowing to calculate a connectivity
probability and connect_feature = feature_mask.
sim_paths = runoutSim(dem = dem,
xy = source_crds,
mu = 0.08,
md = 140,
slp_thresh = 35,
exp_div = 2.5,
per_fct = 1.95,
walks = 1000,
source_connect = TRUE,
connect_feature = feature_mask)
# Plot structure of sim_paths
str(sim_paths)
#> List of 4
#> $ start_cell : num 2803625
#> $ cell_trav_freq: 'table' int [1:2027(1d)] 5 3 14 9 10 1 50 51 33 16 ...
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ all_indices: chr [1:2027] "2774166" "2775633" "2775634" "2775635" ...
#> $ cell_max_vel : num [1:2027] 15.2 12 13.9 16 17 ...
#> $ prob_connect : num 0.859
# Get connectivity probability
sim_paths$prob_connect
#> [1] 0.859
tp_raster <- walksToRaster(sim_paths, dem, method = "cdf_prob")
conn_raster <- connToRaster(sim_paths, dem)
# Crop results
tp_crop <- crop(tp_raster, ext(runout_polygon)+ 500)
conn_crop <- crop(conn_raster, ext(runout_polygon)+ 500)
# Plot results
par(mfrow = c(2,1))
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
mar=c(2,2,1,4), main = "Traverse prob.")
plot(st_geometry(river_channel), col = "lightblue", add = TRUE)
plot(tp_crop, legend = TRUE, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)
plot(crop(hill, ext(runout_polygon)+500), col=grey(0:100/100), legend=FALSE,
mar=c(2,2,1,4), main = "Connectivity prob.")
plot(st_geometry(river_channel), col = "lightblue", add = TRUE)
plot(conn_crop, legend = TRUE, alpha = 0.5, add = TRUE)
plot(st_geometry(runout_polygon), add = TRUE)
plot(source_point, add = TRUE)To simulate runout from multiple source cell locations, we run
runoutSim() iteratively, changing the xy coordinates during
each iteration.
To create a list of xy coordinates from source points (vector), we
use the makeSourceList() function.
# Get coordinates of source points to create a source list object
source_l <- makeSourceList(source_points)
# Perform random walk simulation for each source point
sim_runs <- list()
for(i in 1:length(source_l)){
sim_runs[[i]] <- runoutSim(dem = dem, xy = source_l[[i]],
mu = 0.08,
md = 40,
slp_thresh = 40,
exp_div = 3,
per_fct = 1.9,
walks = 1000)
}We can now convert this list of multiple random walks from different
source cells to raster data using walksToRaster() and
velocityToRaster().
# Coerce results to a raster
trav_freq <- walksToRaster(sim_runs, method = "freq", dem)
vel_ms <- velocityToRaster(sim_runs, dem, method = "max")
trav_prob <- walksToRaster(sim_runs, method = "max_cdf_prob", dem)
# Plot random walks from mulitple source cells
par(mfrow = c(1,3))
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
main = "Traverse frequency")
plot(trav_freq, add = TRUE, alpha = 0.7)
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
main = "Traverse probability (max.)")
plot(trav_prob, add = TRUE, alpha = 0.7)
plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4),
main = "Velocity (m/s)")
plot(vel_ms, col = map.pal("plasma"), add = TRUE, alpha = 0.7)We can also view the results through an interactive leaflet map.
# Create interactive leaflet map
sim_map <-
# Debris flow observations
leafmap(runout_polygons, opacity = 0.3, label = "Runout polygons") %>%
leafmap(source_points, color = '#e74c3c', label = "Source points") %>%
# Terrain data
leafmap(hill, palette = grey(0:100/100), opacity = 1, add_legend = FALSE,
add_image_query = FALSE, label = "Hillshade") %>%
leafmap(dem, palette = viridis::mako(10), opacity = 0.6,
add_image_query = FALSE, label = "DEM") %>%
# Modelling results
leafmap(trav_prob, palette = viridis::viridis(10, direction = -1),
label = "Traverse prob", add_image_query = FALSE) %>%
leafmap(vel_ms, palette = viridis::plasma(10, direction = -1),
label = "Velocity", add_image_query = FALSE)
# Start with mapped zoomed in
sim_map <- leaflet::setView(sim_map, lng = -70.13694, lat = -33.3703, zoom = 13)
# Hide layers
sim_map <- leaflet::hideGroup(sim_map, c("Velocity", "Source points", "DEM", "Hillshade"))
sim_map